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Abstract 

We investigate the dynamical interaction between a galactic disk and surrounding numerous dark 
subhalos as expected for a galaxy-sized halo in the cold dark matter (CDM) models. Our particular 
interest is to what extent accretion events of subhalos into a disk are allowed in light of the observed 
thinness of a disk. Several models of subhalos are considered in terms of their internal density distribution, 
mass function, and spatial and velocity distributions. Based on a series of N-body simulations, we find that 
the disk thickening quantified by the change of its scale height, Azd, depends strongly on the individual 
mass of an interacting subhalo Af su b- This is described by the relation, Azd/Rd — 8^^_ 1 (A/ su b,j /Md) 2 , 
where Rd is a disk scale length, Md is a disk mass, and A" is the total number of accretion events of 
subhalos inside a disk region (< SRd)- Using this relation, we find that an observed thin disk has not ever 
interacted with subhalos with the total mass of more than 15 % disk mass. Also, a less massive disk with 
smaller circular velocity V c is more affected by subhalos than a disk with larger V c , in agreement with the 
observation. Further implications of our results for the origin of a thick disk component are also discussed. 

Key words: cosmology: dark matter — galaxies: formation — galaxies: structure — galaxies: inter- 
actions 



1. Introduction 

The cold dark matter (CDM) paradigm has become a 
standard framework for understanding the structure for- 
mation in the Universe. According to this theoretical 
paradigm, the growing process of self-gravitating struc- 
tures is hierarchical in the sense that small dark matter 
halos virializc first, and aggregate successively into larger 
and larger objects. This clustering process of dark matter 
halos is successful for explaining a wide variety of obser- 
vations including the large-scale distribution of galaxies. 

In this CDM scenario, N-body simulations are an im- 
portant tool in order to investigate the non-linear growth 
of cosmological structures. Early N-body simulations 
based on the CDM models suffered from the so-called 
over-merging problem, i.e., substructures are disrupted 
very quickly within dense environments (Summers et al. 
1995). However, recent high-resolution N-body simula- 
tions have revealed the presence of hundreds of dark mat- 
ter substructures (subhalos) which survive in not only 
cluster scales but also galactic scales (Moore ct al. 1999; 
Klypin ct al. 1999). This large number of subhalos in a 
galaxy-sized halo is in contrast to only about a dozen satel- 
lite galaxies in the Galaxy, which confronts so-called "the 
Missing Satellite Problem" . Several authors have argued 
that this apparent discrepancy could be resolved by con- 
sidering some suppressing process for star formation, such 
as gas heating by an intergalactic ionizing background or 
energy feedback from evolving stars. In whatever models 
relying on the suppression of galaxy formation, a typical 
galaxy-sized halo should contain numerous dark subhalos. 

Then, there is a possibility that a large amount of sub- 



halos interact frequently with a stellar disk embedded in 
the center of a halo, so that the disk would be dynam- 
ically heated and thickened. On the other hand, an ob- 
served galactic disk is rather thin: the scale height (or half 
thickness) is only about ~ 250 pc in the Galaxy. Likewise, 
recent observations of external disk galaxies (Kregel et al. 
2002) suggest that the observed scale height of a disk, 
Zd, is confined to some limiting value relative to the scale 
length of a disk, Rd, i.e, Zd/Rd < 0.2. 

This observed thinness of a disk provides important lim- 
its on the disk heating due to infalling satellites. Toth 
& Ostriker (1992) analytically evaluated this effect and 
concluded that an observed disk like that of the Galaxy 
within the solar radius should have interacted with satel- 
lites with no more than 4 % of the present disk mass 
within the last 5 Gyr. Subsequent numerical simulations 
of an interaction between a disk and a single satellite (e.g. 
Velazquez & White 1999) showed that their analytical 
estimation for the disk heating was somewhat too high 
because an actual interaction process is highly non-linear 
and more complicated than simplified analytical represen- 
tation. Interactions with many subhalos would be much 
more complicated and thus require a more detailed anal- 
ysis. 

Font et al. (2001) have conducted numerical simula- 
tion of interaction between a disk and numerous subhalos 
based on the CDM models. They concluded that the ef- 
fect of subhalos on a disk is rather small, and therefore 
subhalos do not conflict with the presence of a thin disk 
since their orbit seldom take them near the disk. However, 
it is worth noting that in their simulation the initial scale 
height of a disk (700 pc) is already thick compared with 
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the observed one in the Galaxy (~ 250 pc), thereby lead- 
ing to possibly the underestimation of the disk heating 
effect. Their simulation is also limited to only one real- 
ization of subhalos; it is yet unclear whether the derived 
weak effect of subhalos on a disk is general or not. Ardi et 
al. (2003) have investigated more details in this disk heat- 
ing by subhalos. They found that a more massive subhalo 
is more effective to heat the disk than a less massive one. 
However, in their calculation subhalos are represented by 
rigid bodies which never lose their mass irrespective of 
tidal effects of a host galaxy, so that the disk heating is 
overestimated. Also, the applicability of their result to an 
actual disk, especially, to what extent accretion events of 
subhalos into a disk are allowed remains unclear. 

Our aim in this paper is thus to set more useful limits 
on the dynamical interaction between numerous subhalos 
and a galactic disk. For this purpose, we conduct a se- 
ries of numerical simulations, in which a self-gravitating 
disk is embedded in a dark halo containing many subha- 
los. In this work we set an initially thin disk with the 
scale height of 250 pc, in contrast to previous numerical 
studies starting from the scale height of ~ 700 pc much 
larger than the observed one (Velazquez & White f999; 
Font et al. 2001). Several models for the system of sub- 
halos in a host halo are taken into account in terms of 
their mass function, spatial distribution, and velocity dis- 
tribution. We also consider two different models for the 
internal density distribution of subhalos: point-mass and 
extended-mass models. In the latter model, subhalos are 
affected by a tidal field of a host galaxy so that they lose 
their mass in the course of their orbital motions. Based on 
our simulations, we investigate the dependence of the disk 
heating on the model parameters and apply our analysis 
to understanding an observed thin disk in the context of 
the disk heating by subhalos. 

This paper is organized as follows. In § 2 we describe 
our galaxy model which is composed of halo, bulge, and 
disk components. The models of subhalos are also de- 
scribed in this section. In § 3 we present the results of our 
numerical simulations. In § 4 we analyze our results and 
present our prediction for the relation between the disk 
heating by subhalos and an observed thin or thick disk. 
Finally, in § 5 we present our conclusions. 

2. Models 

2.1. The Galaxy Model 

Our galaxy model is composed of three components: a 
disk, a bulge, and a dark halo. To investigate the self- 
gravitating response of the disk component to orbiting 
subhalos, we model the disk by a self-consistent N-body 
realization of stars under the influence of an external force 
provided by the rigid bulge and halo components. The 
methods of Hernquist (1993) are utilized to set up the 
disk consisting of the distribution of N-body particles. 
A detailed description of the technique can be found in 
Hcrnquist's paper. 

The density distribution of the disk is initially axisym- 
metric, pd(R,z), using the cylindrical coordinates (i?, z), 



while the bulge and halo are spherically symmetric, pb{f) 
and Ph{r), respectively, using the galactocentric distance 
r. These density distributions are given by 

p d {R,z) = j—f4z exp(-fl/firf)sech 2 (z/zrf) , (1) 



Pb(r) = 



Ph(r) 



M b Ob 
2tt r(ah + r) 3 ' 
M h a /t exp(-r 2 /r 2 ) 



(2) 



(3) 



2tt 3 / 2 r c r 2 + 7 2 

where M<j, Mb and Mh correspond to the masses of the 
disk, the bulge and the halo, respectively. The disk pa- 
rameters R c i and Zd denote the radial scale length and 
vertical scale height, respectively. The parameter a^ de- 
notes the scale length of the bulge, while 7 and r c are 
the core and cut-off radii for the halo and ah is a nor- 
malization constant. We choose these parameters so that 
the model approximately matches the observed character- 
istics of the Galaxy, and the values of the parameters are 
listed in Tabcl I. It is worth remarking that we consider 
an observed thin scale height for the disk, namely Zd = 245 
pc, in contrast to the models by Font et al. (2001) and 
Velazquez & White (1999) adopting Zd = 700 pc. In or- 
der to prevent the disk from gravitational instability, we 
adopt a stable disk by setting Toomre's Q parameter at 
the solar radius Rq = 8.5 kpc as Q Q = 1.5. The rota- 
tion curve of our galaxy model is shown in Figure 1. The 
rotation speed at the solar radius is V c (Rq) ss 240 km s _1 . 

2.2. Subhalo Models 

We construct a set of subhalo models in our numerical 
simulations, designated as model A to U as tabulated in 
Table 2, to investigate how different physical properties 
of subhalos affect the disk heating process. Each model 
assumption is explained as follows. 

2.2.1. Mass junction, spatial distribution and velocity 
anisotropy 

We consider a mass spectrum for the realization of each 
subhalo with a mass M su b- According to the results of 
cosmological N-body simulations by Moore et al. (1999), 
Klypin et al. (1999), and Ghigna et al. (2000), this mass 
function can be fitted to a power law with an index of 
about —2. We thus adopt the form, 



N(M suh )dM suh cx M~ldM snh 



(4) 



For the convenience of numerical analysis, we set the 
higher and lower mass limits for this mass function desig- 
nated as Mhigh and Mi ow , respectively, and examine the 
role of individual subhalo masses in the disk heating. The 
normalization of equation (4) is given by the total mass 
of the subhalo system, which is about one-tenth of the 
mass of a host halo according to Klypin et al. (1999) and 
Ghigna et al. (2000). We thus set O.lMh as the total mass 
of the subhalo system. 

Recent high-resolution N-body simulations have shown 
that the spatial distribution of subhalos in a host halo is 
less concentrated than the host's density profile (Gao et 
al. 2004), which is often represented by the so-called NFW 
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profile (Navarro, Frenk, & White f997). However, in even 
most recent simulations the mass and force resolutions are 
yet insufficient, so the true spatial distribution of subhalos 
is unclear. In this paper, instead of trying to set a realistic 
spatial distribution (which is yet unknown), we adopt a 
tractable model for it and attempt to extract the general 
results which do not depend on this particular setting. 
Thus, for the initial spatial distribution of subhalos in 
a host halo, we adopt the Hernquist model (Hernquist 
1990), in which the number density n{r) of subhalos at 
the galactocentric distance r is given as 



where a is the scale length in the spatial distribution. It 
is worth noting that the inner density distribution of this 
model is similar to that of the NFW profile. The change 
of the parameter a affects the incidence of subhalo-disk 
interaction, which is mostly effective at r ^ 10 kpc, since 
smaller a yields smaller pericenters and apocenters for the 
orbits of subhalos. This is highlighted in Figure 2, where 
the distributions of the pericenters and apocenters of sub- 
halos are shown for Model F (a = 87.5 kpc) and Model G 
(a = 280 kpc), while having the same velocity distribution 
(sec below). It follows that the number of subhalos orbit- 
ing interior to r 10 kpc is larger for Model F than for 
Model G, and the dependence of this number on the scale 
length a is also seen in other models. 

For the initial velocity distribution of subhalos, we take 
the moments of the collisionless Boltzmann equation fol- 
lowing a procedure described by Hernquist (1993). The 
velocity ellipsoid at each spatial location is calculated from 
the moment equations, and then the velocity components 
are randomly selected from the Gaussian distributions for 
the corresponding velocity ellipsoid. 

We adopt two different models for the velocity 
anisotropy of the subhalo system, which is parameterized 
by j3 = 1 — 0.5(ag+a'^ > )/a^ where ay, a#, and a^ are the ra- 
dial, zenithal, and azimuthal velocity dispersions, respec- 
tively. One is the isotropic model of (3 — 0, which acts as 
our standard model. The other is the radially anisotropic 
model characterized by (3 = 0.5. This anisotropic model is 
motivated by the results of cosmological N-body simula- 
tions (Diemand et al. 2004; Abadi et al. 2006), which show 
the increase of (3 with r, starting (3 ~ at a halo center to 
(3 > 0.5 in its outer parts. For the sake of simplicity, we 
assume (3 is constant along r in our model. 
2.2.2. Effect of Baryon Condensation 

In hierarchical galaxy formation models, stars are 
formed in the condensation of cooled baryon in a halo 
center, subsequently forming a disk component. This con- 
densed baryon or disk pulls the surrounding dark matter 
particles inward, thereby increasing the central concen- 
tration of a dark halo (e.g., Gncdin et al. 2004). This 
effect of baryon condensation is also expected to modify 
the space and velocity distributions of subhalos, compared 
with those obtained by dissipationless N-body simulations 
(Gao et al. 2004). 

We take into account this effect in our model by slowly 



increasing the total masses of the disk and bulge com- 
ponents over a period of 10 Gyr, after setting the initial 
distribution of subhalos in the presence of a (smooth) halo 
alone. When the total masses of the disk and bulge com- 
ponents reach the values listed in Table 1, the position 
and velocity of each subhalo are recorded for the use of 
the further calculations of disk heating. In this experi- 
ment, we treat a subhalo as a point mass and neglect the 
interaction between different subhalos. 

This treatment of the baryon condensation effect is ad- 
mittedly highly ideal and not self-consistent as we neglect 
the simultaneous modification for a smooth halo compo- 
nent 1 . However, the rate of the interactions between sub- 
halos and a disk is somewhat increased by this gravita- 
tional effect of baryonic matter, thereby allowing us to 
carry out a statistically meaningful analysis on the prop- 
erties of the disk heating. In fact, this effect of baryon 
condensation results in a few percent increase in the num- 
ber of subhalos having pericenters smaller than ~ 10 kpc, 
which yields the enough amount of interaction events over 
the interval of numerical simulations. 
2.2.3. Internal Density Distribution 

We consider two different models for the internal den- 
sity distribution of a subhalo: point-mass and extended- 
mass models. In the former models, since point-mass sub- 
halos survive eternally in our simulations, it is postulated 
that subhalos are supplied through their continuous ac- 
cretion into a host halo from outside even if some of them 
disappear due to tidal destruction. In the latter models af- 
fected by tides, we assume a King-model profile character- 
ized by a concentration parameter CKmg = k>g 10 (Rt/R c ), 
where Rt and R c denote tidal and core radii, respectively. 
For these latter models, the tidal effects of the disk (as well 
as the bulge and halo) on subhalos are explicitly taken 
into account. While the adoption of a King-model pro- 
file is admittedly ideal, recent cosmological simulations 
by Kazantzidis et al. (2004) imply that the internal den- 
sity distributions of subhalos may be described reasonably 
well by a more-centrally concentrated universal profile or 
the NFW profile with some tidal outer limit. We therefore 
adjust our King models to match the NFW profile in the 
following manner. Firstly, based on the method outlined 
in NFW, we determine a set of model parameters in the 
NFW profile (see Appendix 1 for details). Secondly, we 
estimate a parameter CKing, whereby the half-mass radius 
of the King model is equal to that of the NFW model. 
Finally, we obtain a tidal radius Rt as a limiting radius of 
the tidal effect of a host galaxy at the initial position of a 
subhalo. Thus, R t is derived from the relation 

M tot (<r) M suh 

—?— = -w (6) 

where Mtot(< r) is the total mass of a host galaxy interior 
to r and Af su b(< Rt) is the mass of a subhalo. For this esti- 
mation of Rt , we assume a spherically symmetric potential 

Our adoption of an isothermal-like profile for a smooth halo 
(equation 3), in comparison with an NFW-likc profile derived 
from N-body simulations, suggests the consideration of baryon 
condensation for the halo setting. 
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Fig. 1. Rotation curve for our disk galaxy model. 

for a host galaxy, where the disk is made spherical with a 
mass distribution Md{r) =Md[l — (l+r/i?d)exp(— r/Rd)]- 
We also take into account the effects of baryon condensa- 
tion for getting the initial position of a subhalo. 

The parameters for point-mass models and extended- 
mass models that wc calculate are summarized in Tabic 
2. 

2.3. Method for Numerical Simulation 

For the point-mass models we use a tree algorithm 
with a tolerance parameter of 9 tor = 0.7 (Barnes & Hut 
1986; Hcrnquist 1987). For the extended-mass models 
we use GRAPE5 systems at the National Astronomical 
Observatory of Japan. The time integration is made with 
the leapfrog method and a fixed time step of 0.41 Myr. 
The softening length for N-body particles is e = 70 pc. 
We use Nd = 46000 particles for the disk and the number 
of the subhalo particles is listed in Table 2. Since the mass 
of the subhalo is negligibly small as compared with that of 
the host galaxy, we neglect the forces between the subhalos 
in the point-mass models. In contrast, for the extended- 
mass models, we fully take into account the gravitational 
interaction between the subhalos for the convenience of 
numerical calculations using GRAPE5. We have followed 
the evolution up to 4.9 Gyr. 

3. Results 

Based on the numerical simulations of the models de- 
fined in the previous section, we examine the effects of 
subhalos on the disk structure and dynamics, especially 
to elucidate the dependence of several different proper- 
ties of subhalos: their internal density distribution, mass 
function, and spatial and velocity distributions. 

In Figure 3 we show the edge-on view of the disk for 
model F at the beginning (t = 0) and the end (t = 4.9 
Gyr) of the simulation. The disk has been thickened and 
tilted by the gravitational interaction with orbiting sub- 
halos. To estimate the change of the disk kinematics and 
thickness at specific locations, wc consider the tilt of the 
disk and use the axes aligned with the principal axes of 
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Fig. 2. Distributions of the pericenter and the apocenter of 
the subhalos in model F (left panel) and model G (right 
panel) . 

the disk inertia tensor. The heating and thickening of the 
disk can be described by the changes of the velocity disper- 
sions (A<7r, Aa z ) and by the increase of the scale height, 
Azrf. To calculate these quantities, the disk is stratified 
in concentric cylindrical annuli with a width of 5R = 700 
pc and the particle properties arc averaged in each annu- 
lus. The scale height in each annulus R of the disk, Zd, 
is defined by the mean square of the z-coordinatcs, i.e., 
Zd(R) =< z 2 > 1 / 2 . We note that at the beginning of the 
calculations, Zd is a constant of 245 pc as given in equation 
(I)- 

In addition to the dynamical effect of subhalos, the sim- 
ulated disk is subject to internal heating due to two-body 
relaxation among the disk particles; this numerical heat- 
ing always takes place in numerical simulations with a 
modest number of particles. We evaluate the effect of in- 
ternal heating by evolving the disk in isolation, i.e., in 
the absence of subhalos. This effect is typically character- 
ized as Azd = 0.23 kpc, Actr = 7.2 km s _1 , and Aa z = 5.7 
km s _1 at the solar radius after 4.9 Gyr. In the followings, 
the notation Azd means the difference of a scale height be- 
tween t = 4.9 Gyr and t = 0, and Actr and Ao~ z also mean 
the change of the velocity dispersions after 4.9 Gyr. 
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3.1. Global Properties of the Disk Evolution 

In Figure 4 we show the kinematical properties of the 
disk for model A, F and G, including the impact of in- 
ternal heating in our calculations. It is evident from this 
figure that thickening of the disk does not occur uniformly 
at all radii; given the complexity of the final disk struc- 
ture, we found it convenient to sample the kinematics at 
R = Rq, 3Rd, and 4i?<j, which is sufficient to provide us 
with a global view of the heating and thickening. The 
growth of the disk thickness in model A is Azd ~ 0.57 
kpc, 0.61 kpc and 0.67 kpc at R — Rq, 3Rd, and ARd, 
respectively, those values in model F are 1.19 kpc, 1.37 
kpc and 1.60 kpc, and those values in model G are 0.28 
kpc, 0.31 kpc, and 0.39 kpc. The increase in the radial 
and vertical velocity dispersions after 4.9 Gyr in model 
A is given as (Aa R ,Acr z ) = (18.8,15.4), (17.4,13.8), and 
(16.6,11.9) km s _1 at R = Rq, 3Rd, and ARd, respec- 
tively, those values in model F are (32.2,26.1), (30.1,23.9), 
and (25.4,22.8) km s" 1 , and those values in model G are 
(9.9,8.2), (10.5,7.5) and (13.5,8.2) km s" 1 . 

In these experiments, the main difference between 
model A and model F resides in individual masses of sub- 
halos parameterized by A/high and M\ ow (see Table 2): for 
model A all subhalos have 10 s M as Mhigh = Mow = 10 8 
Mq, whereas for model F the presence of more massive 
subhalos than 10 s Mq is allowed as Mhigh = 10 9 Mq. 
Therefore the comparison between the results of these two 
models highlights the effect of individual masses of sub- 
halos on the disk heating, where we note that the slight 
difference of the parameter a between the models by a fac- 
tor 1.25 yields essentially no difference in the results. It 
is clear that the presence of a few, but massive subhalos 
(model F) is more effective for the disk heating than the 
case of many but less massive ones (model A), as already 
pointed out by Ardi et al. (2003). This suggests that the 
disk heating process is more sensitively enhanced than 
being proportional to individual subhalo masses; massive 
subhalos are more important for the disk heating. Also, 
in comparison with model F, model G with the same val- 
ues for Mhigh and M\ ow yields a weak effect on the disk. 
The main difference between these two models is the spa- 
tial distribution of subhalos parameterized by a (a = 87.5 
kpc and 280 kpc for model F and G, respectively), which 
affects the pericenter distributions of subhalos especially 
at r < 10 kpc (see Figure 2). Therefore, we find that the 
number of subhalos which cross the disk is also important 
in quantifying the disk heating. 

Figure 5 shows the growth of disk velocity dispersions 
in the radial and vertical directions at R = Rq for model 
F (the point-mass model), K and L (the extended-mass 
model). This figure shows that the disk velocity dispersion 
for model F continues to grow throughout the simulations, 
whereas for model K and L the growth of the velocity 
dispersion almost stops at t ~ 1 Gyr. It is worth noting 
that in these latter models subhalos can lose their mass 
at the interaction with the disk, unlike the former point- 
mass models in which subhalo masses remain the same. 
Thus, for the extended-mass models the effect of subhalos 
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Fig. 4. The kinematical properties of the disk after 4.9 Gyr 
for different subhalo models. Thick solid lines show the initial 
setting and thin solid lines show the effects of internal heating 
in the absence of subhalos. Dashed, dotted, and dash-dotted 
lines show model A, F, and G, respectively. 
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Fig. 3. Growth of the disk thickness for model F. The left and right panels show the egde-on view of the disk at the beginning 
(t = 0) and the end (t = 4.9 Gyr) of the simulation, respectively. 



on the disk is temporal; subhalos appear to lose almost 
all of their mass at the first interaction with the disk, so 
that their role in the disk heating is effective only at the 
first interaction with the disk. 

3.2. Disk Thickness vs. Subhalo Masses 

In this section we investigate in more detail the effect 
of orbiting subhalos on the growth of the disk thickness. 
As discussed above, the change of the disk scale height, 
Zd, depends on R, in such a manner that Azd is some- 
what larger at larger R. To quantify this disk thickening 
as a function of accreted subhalo masses, we utilize the 
disk scale height at the outer edge of the disk, R = R ou t, 
where the stellar light distribution is expected to diminish 
steeply with R. This truncation of a stellar disk has actu- 
ally been observed, where the stellar light declines more 
steeply than an exponential profile for a main disk and 
drops to low values beyond the so-called truncation ra- 
dius (van der Kruit & Searle 1981a, b, 1982; Kregel et 
al. 2002). This usually occurs at a radius of 3 — 5 disk 
scale length (Kregel ct al. 2002), so in our work we adopt 
^?out = 3i?d as a characteristic outer edge of the disk, at 
which the change of the disk scale height is evaluated. 
Also, to analyze the relation between the change of the 
disk thickness and the orbits of the accreted subhalos, we 
estimate the number of times which each subhalo crosses 
the disk region at R < 3Rd in the course of its orbital 
motion. 

We have carried out the simulations for several different 
models of point-mass subhalos with j3 = 0. Figure 6 shows 
the relation between the growth of the disk scale height 
at R = 3Rd (i.e. Azd/Rd) and the combination of each 
mass of a subhalo M Bn b,i and the number of times which 
it crosses the disk at R < 3Rd in the course of its orbital 
motion, denoted as Ni [i.e. '}2 i Ni(M su \ 1 ,i/Md) for the left 
panel and ^2 i N i (M su ^ i /Md) 2 for the right panel]. Here 
the abscissa for the left panel is proportional to the to- 
tal accreted mass of subhalos, whereas that for the right 
panel is proportional to the sum of the squared masses 
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Fig. 5. Growth of the disk velocity dispersion in the radial 
and vertical directions at R = Rq for model F (dashed lines), 
K (dotted lines), and L (dash-dotted lines). In model F sub- 
halos are represented by point masses, whereas in model K 
and L they have a King-model profile. Thick solid lines show 
the effect of internal heating in the absence of subhalos. 
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of subhalos. As is evident, while the left panel shows no 
correlation between the abscissa and ordinate axes, the 
right panel shows a significantly tight correlation, thereby 
indicating that the growth of the disk thickness is pro- 
portional to the sum of the squared masses of subhalos; 
if so, the disk thickening is more enhanced for more mas- 
sive individual subhalos as already shown in the previous 
subsection. We thus investigate this relation for different 
subhalo models as shown in Figure 7. In the left panel 
we consider the extended-mass models with j3 = (filled 
squares) in comparison with the point-mass models with 
(3 — (open circles) . We note here that in the extended- 
mass models subhalos lose almost all of their mass at the 
first interaction with the disk, so we adopt N. \ = 1 for such 
a case. The right panel shows the case of an anisotropic 
velocity distribution with (3 = 0.5 for the point-mass mod- 
els (filled triangles) and with (3 = for the point-mass 
models (open circles). As is evident from these panels, 
several different models yield an almost universal relation 
between Az d /R d at R = 3R d and £\ Ni(M suh ^/M d ) 2 at 
R < 3R d . 
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Fig. 7. The same as Figure 6 but for several different subhalo 
models plotted in logarithmic scales. The left panel shows the 
extended-mass models with /3 = (filled squares) in compar- 
ison with the point- mass models with (3 = (open circles), 
whereas the right panel shows the case of an anisotropic ve- 
locity distribution with (3 = 0.5 for the point-mass models 
(filled triangles) and point-mass models with fi = (open cir- 
cles). The solid lines show the fitting to the relations for the 
point-mass models with /3 = by the least-squares method, 
dashed lines for the extended-mass models with /3 = 0, and 
dotted lines for the point-mass models with j3 = 0.5. 



4. Discussion 

4-1- Dependence of the Disk Heating on Subhalo Masses 

From the results of § 3.2, we find that the disk thickness 
is increasing with the accretion of subhalos into a disk. 
Our numerical experiments suggest the following universal 
relation, 



^ N, (M| / M d r 



Az d 

Rd 



AT ( -^SUb,j 



(7) 



Fig. 6. Relation between the growth of the disk scale height 
at R = 3Rd and the combination of each mass of a subhalo 
M su i, i and the number of times which it crosses the disk at 
R < 3Rd, denoted as JVj, for the point- mass subhalo models 
with p = (model A to J listed in Table 2). The abscis- 
sas in the left and right panels show y\ Nj(M su ^ i/M l i) and 
. -^(Maub.j/Mrf) 2 , respectively. Notice that the effect of 
internal heating has been subtracted in these plots. 



where a is a constant of ~ 8, R d is the disk scale length, 
z d is the disk scale height at R = 3R d , and iVj is the 
number of times that subhalos with an individual mass 
of M su b,i cross a disk at R < 3R d (noting that JVi = 1 for 
the extended-mass models as subhalos lose their mass at 
their first interaction with the disk). In this expression, 
a subscript i denotes an individual subhalo given at t = 
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in our simulation. An alternative, more useful expression 
based on equation (7) is derived as follows. Supposing 
that the disk has experienced the accretion of subhalos 
having an individual mass of Af su b.j with j = 1,...,N at 
R < 3i?d, where the repeated accretion of a subhalo in the 
course of its orbital motion is regarded as a separate ac- 
cretion event with a mass M su b, j , and N denotes the total 
number of such events. Then, we obtain the relation, 



R d 



\ ^ / A/ SU b,j 
M d 



(8) 



which holds a more useful form for any applications than 
equation (7). Notice that although this relation is derived 
from the simulations over the interval of 4.9 Gyr, this 
is applicable to longer time evolution by considering the 
total number of subhalos crossing a disk, N. 

As equation (8) indicates, the increase of the disk thick- 
ness is proportional to the square of the masses of accreted 
subhalos. This mass dependency in the disk heating pro- 
cess can be understood if we consider the transfer of ki- 
netic energy from the subhalos to the disk through dy- 
namical friction. Here we present the summary of this 
derivation and more details are shown in Appendix 2. 

Firstly, the vertical equilibrium between kinetic and 
gravitational energy allows us to relate the velocity dis- 
persion o z and the scale height Zd of a disk. Assuming 
that a disk is an isothermal sheet, we obtain 



ct 2 = 2TTGY,(R)z d , 



(9) 



where S(i?) is a surface density of a disk at R (e.g. Spitzer 
1942). Secondly, we consider the energy input into disk 
stars getting through subhalo-disk interaction: the energy 
loss of a subhalo is equal to the energy pumped into a disk. 
This energy loss AE SU ^ is derived by the integral over an 
orbit of a subhalo, 



AK 



(10) 



where -Fdrag corresponds to a dynamical friction. Using 
the Chandrasckhar formula for -Fdrag, each subhalo with 
a mass M su b is subject to a frictional force with i*d ra g °c 
M 2 ub . Finally, combining equation (9) with (10), we ob- 
tain, 



Az d oc Act? oc AE s . 



Fdra 



Zd oc Af 2 ub 



(11) 



Therefore, the dependence of Az d on M su b as we have 
obtained in § 3 is understood in the framework of a dy- 
namical friction between a subhalo and a disk. 

4-2. Comparison with an Observed Thin Disk 

Recent observations of external disk galaxies by Kregel 
ct al. (2002) have suggested that the thickness of a (thin) 
disk is confined to some limiting value relative to a scale 
length of a disk, which is expressed as Zd/Rd < 0.2. Kregel 
et al. (2002) have also shown that the distribution of 
Zd/Rd tends to have an increasing dispersion with increas- 
ing maximum circular velocity V c of a disk, in such a man- 
ner that larger V c allows smaller z d / Rd\ conversely, a disk 
with smaller V c is likely thicker. 



An observed thin disk with Zd/Rd < 0.2 suggests that 
the accretion of subhalos has been rather insignificant 
since a disk with a current mass was formed. Using 
equation (8), this observed limit implies (J2j M^ j) 1 / 2 < 
Q.15M d - Thus, we find that an observed thin disk has not 
ever interacted with subhalos with the total mass of more 
than 15 % disk mass. 

The dependence of Zd/Rd on V c may be understood 
as follows. Let (M 2 uh ) as the mean square of a subhalo 

mass, defined as <M 2 ub ) = £f M 2 uhj /N. Then, equation 
(8) can be written as 



Rd M 2 d 



(12) 



We suppose that N, the number of subhalos which cross a 
disk at R < 3R d , is roughly proportional to the spherical 
volume with a radius 3Rd inside a dark halo, given the spa- 
tial distribution of subhalos. This reads TV oc R d - Also sup- 
pose that the central disk surface density S = M d /{2nR d ) 
is nearly the same for a disk galaxy, which may correspond 
to nearly the same central surface brightness for a bright 
disk (e.g., Freeman 1970). Then equation (12) is rewritten 
as, 



^ oc <M 2 ub ) • V- 



(13) 



where V c is a disk circular velocity derived from V c oc 
GMd/Rd- 

Thus, if the observed disk thickness is controlled by 
disk-subhalo interaction, i.e., Zd ~ Azd, and (Ajf 2 ub ) is 
roughly the same for each disk galaxy, then we find that 
the effect of subhalos on a disk with a smaller mass is 
more significant than a disk with a larger mass. This is in 
good agreement with the observation (Kregel et al. 2002) 
that a disk with smaller V c is likely thicker. 

4-3. The Relation to the Origin of a Thick Disk 

In recent years, new datasets for the thick disk com- 
ponent of the Galaxy as well as for a thick disk of an 
external galaxy have been available, showing several im- 
portant properties of a thick disk. Based on the third 
data release of the Sloan Digital Sky Survey (York et al. 
2000), Allcnde Prieto et al. (2006) have found that the 
stars belonging to the thick disk have no vertical metal- 
licity gradient. The rotational velocity for the same stars 
however shows a vertical gradient of ~ —16 km s _1 kpc -1 
between 1 and 3 kpc from the Galactic plane. The obser- 
vations of an external disk galaxy (Yoachim & Dalcanton 
2006) have shown that the ratio of the total luminosity of 
a thick disk to that of a thin one is related to the circular 
velocity of a galaxy, in such a manner that the ratio tends 
to increase with decreasing circular velocity; a less mas- 
sive galaxy has a brighter (and possibly more massive) 
thick disk relative to a thin one. These properties of a 
thick disk are expected to set important constraints on its 
origin. 

Here, we consider a possibility that a thick disk has been 
formed by the dynamical effect of numerous dark subha- 
los on a pre-existing disk. This interaction effect may be 
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rather weak at the current epoch as only a small fraction 
of subhalos would cross a disk (Font et al. 2001), but the 
effect at early times may have been more significant due to 
a smaller disk mass and larger accretion rate of subhalos 
onto a parent halo. If so, a pre-existing disk should have 
suffered from considerable heating at early times, which 
results in the formation of a thick disk; subsequent slow 
accumulation of baryonic gas in a plane may form a thin 
disk component. 

To assess this possibility in light of the observed prop- 
erties of a thick disk, we investigate the properties of the 
disk thickened by subhalos in our model. Figure 8 shows 
the distribution of the stars in model F at the end of the 
simulation (t = 4.9 Gyr), differentiated by the initial (t = 0) 
positions at either z < 245 pc (left panel) or z > 245 pc 
(right panel). As is evident, the disk stars are well mixed 
by the disk heating. This suggests that even if a pre- 
existing disk had a metallicity gradient, the disk heating 
by interactions with subhalos would have wiped out this 
gradient, leaving a thick disk in agreement with the obser- 
vations. It is also suggested that this disk heating process 
prompts a vertical gradient in rotational velocity. This is 
explained as follows. A heated disk puffs up not only in 
the vertical but also in the radial direction, where the stars 
located at high \z\ have large radial velocity dispersion 
compared with those at low \z\ and thus have small rota- 
tional velocities on average due to the effect of asymmetric 
drift: given a gravitational force inward, the increase of 
radial pressure of stars reduces the effect of a centrifugal 
force. This vertical gradient in rotational velocity is actu- 
ally obtained in our simulation models as shown in Figure 
9. It is found that the gradient amounts to — 10 ~ —30 
km s _1 kpc" 1 between 1 and 3 kpc from the plane, in 
good agreement with the observations. This experiment 
suggests that the presence of a vertical gradient in rota- 
tional velocity of a thick disk may be an important clue 
to distinguishing the scenarios for the origin of a thick 
disk; models invoking monolithic disk collapse (Burkert, 
Truran & Henslcr 1992) or chaotic merging event of build- 
ing blocks (Brook et al. 2004) at early times of a galaxy 
may have difficulties in this regard. 

More definite picture for the formation of a thick disk 
must await more elaborate modeling of a forming galaxy 
in the context of hierarchical clustering. In particular, 
each galaxy has a different merging history of subhalos 
and so different spatial and velocity distributions, which 
inevitably affects the interpretation for the observed prop- 
erties of a thick disk, such as for its fraction as a function 
of a disk mass. Also, the observations show that almost 
all of disk galaxies have a distinct thick disk in addition 
to a thin component. This implies that the dynamical 
effects of subhalos on a disk were significant only before 
a specific epoch and the subsequent formation of a thin 
disk has been unaffected by subhalos, although it is yet 
unclear if it is applicable to all galaxy-sized halos having 
different merging histories. Therefore, to assess the sce- 
nario, more detailed numerical studies are required, tak- 
ing into account the growth of both a dark halo and disk 
simultaneously. 
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Fig. 8. Distribution of stars in model F at the end of the 
simulation (t = 4.9 Gyr), differentiated by their initial (t = 0) 
positions at z < 245 pc (left panel) and z > 245 pc (right 
panel). 
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5. Conclusions 

We summarize our conclusions as follows: 

• The dynamical effects of subhalos on a disk are rep- 
resented by the relation between the change of the 
disk scale height Az d (measured at the disk edge 
R = 3R d ) and individual masses of subhalos M su b, 
i.e., Az d /R d ~ 8Y,f=i( M snh,j/M d ) 2 , where R d is a 
disk scale length, M d is a disk mass, and N is the 
total number of accretion events of subhalos inside 
a disk region (< 3R d ). 

• If subhalos with the total mass of more than 15 % 
disk mass interact with a disk, then the disk thick- 
ness is made larger than the observed range. 

• A less massive disk with smaller circular velocity 
V c is found to be more affected by subhalos than a 
disk with larger V c , which is in agreement with the 
observed properties of a thin disk. 

• Stars in a significantly thickened disk by subhalos 
appear to be well mixed and show a vertical gradi- 
ent in their rotation velocity, being similar to the 
observed properties of the thick disk in the Galaxy. 

We note that the relation (8) we have obtained here is 
universal and thus useful for the applications to any rele- 
vant issues, including the dynamics of an evolving stellar 
disk at the center of a growing dark halo. Such detailed 
studies of a galactic disk in comparison with recently in- 
creasing datasets of a remote disk galaxy will be of great 
importance and is left to future work. 

The numerical computations reported here were car- 
ried out on GRAPE systems (project ID: g05b05) kindly 
made available by the Astronomical Data Analysis Center 
(ADAC) at the National Astronomical Observatory of 
Japan (NAOJ). This work has been supported in part by a 
Grant-in- Aid for Scientific Research (15540241, 17540210) 
from the Japanese Ministry of Education, Culture, Sports, 
Science and Technology. 
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Appendix 1. NFW profile 

The NFW density profile is given by 

P = Petit i , s M . TT2 ' ( A1 ) 

(r/r s )(l + r/r s y 

where p c ; r t is the critical density of the Universe, r s is a 
scale radius, and Sq is a characteristic density contrast. 
Following NFW, we define the limiting radius of a virial- 
ized halo, r2oo, to be the radius within which the mean 
mass density is 200p cr i t . Also, the concentration param- 
eter of a halo is defined as c = r2oo / r a , with which <5o is 
given as 



200 



3 [ln(l + c)-c/(l + c)] 



(A2) 



To put the NFW density profile in a cosmological con- 
text, we need to calculate the concentration factor c, 
which is related to 5q via equation (A2). The appropri- 
ate value of c depends on halo formation history and on 
cosmology. NFW proposed a simple model for c based on 
halo formation time. The formation redshift z co n of a halo 
identified at z — with mass M is defined as the redshift 
by which half of its mass is in progenitors with mass ex- 
ceeding fM, where / is a constant. With this definition, 
z co n can be computed by simply using the Press-Schechter 
formalism (e.g. Lacey & Cole 1993), 



erfc ■ 



S c (z con )-S° c 



[y/2[A*(fM)-Al(M)] 



(A3) 



where Aq(M) is the linear variance of the power spectrum 
at z = smoothed with a top- hat filter of mass M, 8 c (z) 
is the density threshold for spherical collapse by redshift 
z, and 6® = 5 C (Q). NFW found that the characteristic 
overdensity of a halo at z = is related to its formation 
redshift z co \\ by 



No. 
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S (M, /) = C(/)O [l + 2coii(M, /)] s 



(A4) 



where the normalization C (/) depends on / and ilo is the 
current density parameter of the Universe. We will take 
/ = 0.01 as suggested by the N-body results of NFW. In 
this case C(/) «3x 10 3 . Thus, for a halo of given mass 
at z = 0, one can obtain the concentration factor c from 
equations (A2)-(A4). In practice, we first solve z co n from 
equation (A3) and insert the value of z co \\ into equation 
(A4) to get 8q. We then use this value of i5o in equation 
(A2) to solve for c. 

In this experiment, we adopt a standard set of cosmolog- 
ical parameters as Sl = 0.3, A = 0.7, h = 0.7, and eg = 1.3, 
where A is a cosmological constant, h is a normalized 
Hubble constant of h = i/o/100 km s -1 Mpc -1 , and ag 
parameterizes density fluctuations at 8 hr 1 Mpc. 

Appendix 2. Derivation of equation (8) 

We show here the derivation of equation (8) based on 
the following dimensional analysis. 

Firstly, we derive the relation between the vertical ve- 
locity dispersion <j z of disk stars and the scale height Zd 
of a self-gravitating axisymmetric disk in cylindrical co- 
ordinates (R,z). For the limit of a thin disk where z- 
derivatives dominate over i?-derivatives, an equation for 
vertical equilibrium and a Poisson equation read, respec- 
tively, 

!h 2 



1 d(pd<T 2 z 
dz 



Pd 



oz 



■53- - 4nG P > 



(A5) 



where pd is a mass density of a disk and $ is a gravitational 
potential. Integrating the second equation over all z gives 
the surface mass density vs. $, i.e., 2ttCS(R) ~ 

d^/dz. Then, inserting this into the first equation, and 
assuming pd oc cxp(— z/zd) and a z is constant, we obtain 
Zd — cr 2 /2ttGT,(R) . This equation suggests that the change 
of the disk thickness is related to the change of the velocity 
dispersion, i.e., 



Az d 



Ac 



(A6) 



2ttGE(R) ' 

Secondly, the change of the velocity dispersion Act 2 is 
related to the energy input into disk stars getting through 
subhalo-disk interaction: the energy loss of a subhalo is 
equal to the energy pumped into a disk. Denoting this 
energy loss as AE sn b per a subhalo, the resultant Act 2 for 
a disk with a total mass Md reads, 

Act 2 = ^ . (A7) 

Md 

We note that Ai? su b is derived by the integral over an 
orbit of a subhalo, 



(A8) 



where Fkrag corresponds to a dynamical friction. Using 
the Chandrasckhar formula for Fdrag, each subhalo with 
a mass M su b is subject to a frictional force with 



^drag — 



47rG 2 M s 2 ub p d lnA 



erf (X) - ^Le- x2 

\/7T 
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(A9) 



where In A is the Coulomb logarithm, X = v su b/\/2<j with 
a being the disk velocity dispersion, and v su h is the sub- 
halo's velocity. For u su b, we suppose that it is represented 
by a virial velocity of a smooth dark halo with a mass Mh 
and a radius rjj, giving w 2 ub ~ GMh /th ■ Furthermore, if 
Md and Rd for a disk is some fraction of Mh and rjj, 
given as Md = fiMjj and Rd = fiTH, where typically 
fi ~ h ~ 0(10-!), we obtain v 2 suh ~ GM d /Rd- For p d , 
we set p d ~ M d /Rjz d . 

Finally, using the above equations, the change of the 
disk thickness induced by accretion events of subhalos 
is estimated by, 



NA( ?z 

NM 2 u 

sub 



NG 2 z d M 2 uh p d /v 2 



sub 



GEM d 



oc 



XM d R d 
Thus, we obtain 

which is consistent with equation (8). 



(A10) 



(All) 
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Table 1. Galactic parameters 





Symbol 


Value 


Disk: 








N d * 


46000 




M d 


5.6 x 10 10 M Q 




Rd 


3.5 kpc 




Zd 


245 pc 




Qe 


1.5 




Rq 


8.5 kpc 




e 


70 pc 


Bulge: 








M b 


1.87 x 10 10 M 




Cbb 


525 pc 


Halo: 








Mh 


7.84 x 10 11 M 




7 


3.5 kpc 




r c 


84 kpc 



* The number of particles used for the disk. 



Table 2. The parameters of the subhalos 



Model 



Number of 
subhalos 



a 
[kpc] 



Mhigi 
[M Q ] 



Mlow 

[Mq] 



^sub 



point-mass models with (3 = 



A 


784 


70 


10" 


10 8 




B 


784 


140 


10 8 


10 8 




C 


392 


140 


2 x 10 8 


2 x 10 8 




D 


261 


140 


3 x 10 8 


3 x 10 8 




E 


200 


175 


4x 10 8 


4x 10 8 




F 


318 


87.5 


10 9 


10 8 




G 


313 


280 


10 9 


10 8 




H 


175 


140 


10 10 


10 s 




I 


1141 


175 


10 10 


10 7 




J 


1959 


140 


10 9 


10 7 






extended- 


■mass mo 


dels with (3 = 




K 


318 


87.5 


10 y 


10 8 


182 


L 


175 


140 


10 10 


10 8 


182 


M 


362 


24.5 


10 9 


10 8 


182 


N 


200 


175 


4x 10 8 


4x 10 8 


182 





112 


70 


7x 10 8 


7x 10 8 


182 


Pt 


280 


52.5 


10 10 


10 8 


170 


Q 


173 


24.5 


10 10 


10 8 


170 


point-mass models with (3 = 0.5 


R 


197 


140 


10 1U 


10* 




S 


362 


87.5 


10 9 


10 8 




T 


249 


140 


3 x 10 9 


10 8 




U 


361 


157.5 


10 9 


10 8 





* The number of particles used tor each subhalo. 

t Only in this model, the total mass of the subhalo system is 13 % 

of that of the host galaxy. 



